Points
snxy <- data.frame(name = "Space Needle", x = -122.3493, y = 47.6205)
# x:longitude, y:latitude
space_needle <- st_as_sf(snxy, coords = c("x", "y"), crs = 4326)
print(space_needle)
## Simple feature collection with 1 feature and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: -122.3493 ymin: 47.6205 xmax: -122.3493 ymax: 47.6205
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## name geometry
## 1 Space Needle POINT (-122.3493 47.6205)
st_coordinates(space_needle)
## X Y
## 1 -122.3493 47.6205
shxy <- data.frame(name = "Savery Hall", x = -122.3083, y = 47.6572)
savery_hall <- st_as_sf(shxy, coords = c("x", "y"), crs = 4326)
# rbind() to put two points in one data frame
pts <- rbind(space_needle, savery_hall)
print(pts)
## Simple feature collection with 2 features and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: -122.3493 ymin: 47.6205 xmax: -122.3083 ymax: 47.6572
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## name geometry
## 1 Space Needle POINT (-122.3493 47.6205)
## 2 Savery Hall POINT (-122.3083 47.6572)
plot(pts$geometry, axes = TRUE)

Lines
# create a linestring sf data frame
lnstr <- st_sfc(st_linestring(st_coordinates(pts)), crs = 4326)
lnstr <- as_tibble(lnstr) %>% mutate(od = "Space Needle, Savery Hall")
plot(pts$geometry, axes = TRUE)
text(x = st_coordinates(pts), labels = pts$name)
plot(lnstr$geometry, col = 2, add = TRUE)
## Polygons
zooxy <- data.frame(name = "Woodland Park Zoo", x = -122.3543, y = 47.6685)
wp_zoo <- st_as_sf(zooxy, coords = c("x", "y"), crs = 4326)
# rbind() to put two points in one data frame
pts <- rbind(pts, wp_zoo)
(plygn <- st_sfc(st_polygon(list(st_coordinates(rbind(pts, space_needle)))), crs = 4326))
## Geometry set for 1 feature
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -122.3543 ymin: 47.6205 xmax: -122.3083 ymax: 47.6685
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## POLYGON ((-122.3493 47.6205, -122.3083 47.6572,...
plot(plygn, col = "cyan", axes = TRUE)
plot(lnstr$geometry, col = 2, add = TRUE, lwd = 3)
plot(pts$geometry, add = TRUE, cex = 2)
text(x = st_coordinates(pts), labels = pts$name)
## Importing spatial data sets
# path to the data
mydatadir <- file.path("C:", "users", Sys.getenv("USERNAME"), "Documents","study", "UW", "GIS_R","data")
zippolyfname <- file.path(mydatadir, "zip_poly.gdb")
# avoid reading over and over
if(!exists("zipcodes")){
zipcodes <- st_read(dsn = zippolyfname, layer = "zip_poly", as_tibble = TRUE, geometry_column = "Shape")
}
## Reading layer `zip_poly' from data source `C:\Users\Yohan_Min\Documents\study\UW\GIS_R\data\zip_poly.gdb' using driver `OpenFileGDB'
## Simple feature collection with 30924 features and 8 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -179.1473 ymin: 13.23419 xmax: 179.7785 ymax: 71.39048
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
# change the data frame's column names to lowercase
colnames(zipcodes) <- tolower(colnames(zipcodes))
# after renaming columns it is necessary to re-establish which column contains the geometry
st_geometry(zipcodes) <- "shape"
zip_wa <- zipcodes %>% filter(state == "WA")
head(zip_wa)
## Simple feature collection with 6 features and 8 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -124.0031 ymin: 46.0401 xmax: -119.2583 ymax: 49.0014
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 9
## zip_code po_name state population pop_sqmi sqmi shape_length shape_area
## <fct> <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 00072 Nation~ WA -99 -99 1.22e3 4.20 0.386
## 2 00073 Usdoe ~ WA -99 -99 3.24e2 1.63 0.0984
## 3 00074 Yakima~ WA -99 -99 1.35e3 4.51 0.408
## 4 00076 Okanog~ WA -99 -99 7.31e2 3.83 0.232
## 5 00195 Long I~ WA -99 -99 8.37e0 0.536 0.00254
## 6 98001 Auburn WA 35721 2003. 1.78e1 0.628 0.00549
## # ... with 1 more variable: shape <MULTIPOLYGON [°]>
plot(x = zip_wa$shape, axes = TRUE)

hospitals <- st_read(file.path(mydatadir, "medical_facilities/medical_facilities.shp"))
## Reading layer `medical_facilities' from data source `C:\Users\Yohan_Min\Documents\study\UW\GIS_R\data\medical_facilities\medical_facilities.shp' using driver `ESRI Shapefile'
## Simple feature collection with 154 features and 14 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 1254366 ymin: 78085.2 xmax: 1395044 ymax: 287325.1
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=47.5 +lat_2=48.73333333333333 +lat_0=47 +lon_0=-120.8333333333333 +x_0=500000.0000000001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
h2o <- st_read(file.path(mydatadir, "wtrbdy/wtrbdy.shp"))
## Reading layer `wtrbdy' from data source `C:\Users\Yohan_Min\Documents\study\UW\GIS_R\data\wtrbdy\wtrbdy.shp' using driver `ESRI Shapefile'
## Simple feature collection with 15838 features and 16 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 607316.1 ymin: -256676.3 xmax: 1617446 ymax: 765087.4
## epsg (SRID): NA
## proj4string: +proj=lcc +lat_1=47.5 +lat_2=48.73333333333333 +lat_0=47 +lon_0=-120.8333333333333 +x_0=500000.0000000001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
# type = "n" not to plot, but sets xlim and ylim
plot(hospitals$geometry, type = "n", axes = TRUE)
# water
plot(h2o$geometry, col = "cyan", border = 0, add = TRUE)
# ZIP code areas
plot(zip_wa$shape, add = TRUE, col = 0, border = 1)
# hospital points
plot(hospitals$geometry, add = TRUE)
box()

# different coordinates
# st_crs(hospitals)
# st_crs(zip_wa)
Exporting data
# st_write(obj = zip_wa, dsn = file.path(mydatadir, "zip_wa.shp"))
st_write(obj = pts, dsn = file.path(mydatadir, "r_gis.gpkg"), layer = "pts")
## Updating layer `pts' to data source `C:/users/Yohan_Min/Documents/study/UW/GIS_R/data/r_gis.gpkg' using driver `GPKG'
## Updating existing layer pts
## Writing 3 features with 1 fields and geometry type Point.
st_write(obj = lnstr, dsn = file.path(mydatadir, "r_gis.gpkg"), layer = "lnstr")
## Updating layer `lnstr' to data source `C:/users/Yohan_Min/Documents/study/UW/GIS_R/data/r_gis.gpkg' using driver `GPKG'
## Updating existing layer lnstr
## Writing 1 features with 1 fields and geometry type Line String.
st_write(obj = plygn, dsn = file.path(mydatadir, "r_gis.gpkg"), layer = "plygn")
## Updating layer `plygn' to data source `C:/users/Yohan_Min/Documents/study/UW/GIS_R/data/r_gis.gpkg' using driver `GPKG'
## Updating existing layer plygn
## Writing 1 features with 0 fields and geometry type Polygon.
st_write(obj = zip_wa, dsn = file.path(mydatadir, "r_gis.gpkg"), layer = "zip_wa")
## Updating layer `zip_wa' to data source `C:/users/Yohan_Min/Documents/study/UW/GIS_R/data/r_gis.gpkg' using driver `GPKG'
## Updating existing layer zip_wa
## Writing 547 features with 8 fields and geometry type Multi Polygon.
Coordinates
snxy <- data.frame(name = "Space Needle", x = -122.3493, y = 47.6205)
space_needle <- st_as_sf(snxy, coords = c("x", "y"))
st_crs(space_needle)
## Coordinate Reference System: NA
st_crs(space_needle) <- 4326
st_crs(space_needle)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.4-4, (SVN revision 833)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/Yohan_Min/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/Yohan_Min/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.3-1
epsg <- make_EPSG()
utm10 <- epsg[grep("UTM.*10", epsg$note),]
kable(utm10) %>% kable_styling(bootstrap_options = "striped", full_width = F, position = "left")
|
|
code
|
note
|
prj4
|
|
1647
|
3157
|
# NAD83(CSRS) / UTM zone 10N
|
+proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
|
|
2207
|
3717
|
# NAD83(NSRS2007) / UTM zone 10N
|
+proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
|
|
2230
|
3740
|
# NAD83(HARN) / UTM zone 10N
|
+proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
|
|
3031
|
6339
|
# NAD83(2011) / UTM zone 10N
|
+proj=utm +zone=10 +ellps=GRS80 +units=m +no_defs
|
|
4126
|
26710
|
# NAD27 / UTM zone 10N
|
+proj=utm +zone=10 +datum=NAD27 +units=m +no_defs
|
|
4274
|
26910
|
# NAD83 / UTM zone 10N
|
+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs
|
|
4850
|
32210
|
# WGS 72 / UTM zone 10N
|
+proj=utm +zone=10 +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.2263 +units=m +no_defs
|
|
4910
|
32310
|
# WGS 72 / UTM zone 10S
|
+proj=utm +zone=10 +south +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.2263 +units=m +no_defs
|
|
4970
|
32410
|
# WGS 72BE / UTM zone 10N
|
+proj=utm +zone=10 +ellps=WGS72 +towgs84=0,0,1.9,0,0,0.814,-0.38 +units=m +no_defs
|
|
5030
|
32510
|
# WGS 72BE / UTM zone 10S
|
+proj=utm +zone=10 +south +ellps=WGS72 +towgs84=0,0,1.9,0,0,0.814,-0.38 +units=m +no_defs
|
|
5091
|
32610
|
# WGS 84 / UTM zone 10N
|
+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
|
|
5159
|
32710
|
# WGS 84 / UTM zone 10S
|
+proj=utm +zone=10 +south +datum=WGS84 +units=m +no_defs
|
|
5470
|
6653
|
# NAD83(CSRS) / UTM zone 10N + CGVD2013 height
|
+proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +vunits=m +no_defs
|
kable(epsg %>%filter(code == 2927)) %>% kable_styling(bootstrap_options = "striped", full_width = F, position = "left")
|
code
|
note
|
prj4
|
|
2927
|
# NAD83(HARN) / Washington South (ftUS)
|
+proj=lcc +lat_1=47.33333333333334 +lat_2=45.83333333333334 +lat_0=45.33333333333334 +lon_0=-120.5 +x_0=500000.0001016001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
|
(space_needle_utm10 <- space_needle %>% st_transform(26910))
## Simple feature collection with 1 feature and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: 548894.1 ymin: 5274327 xmax: 548894.1 ymax: 5274327
## epsg (SRID): 26910
## proj4string: +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## name geometry
## 1 Space Needle POINT (548894.1 5274327)
(space_needle <- space_needle %>% st_transform(26910))
## Simple feature collection with 1 feature and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: 548894.1 ymin: 5274327 xmax: 548894.1 ymax: 5274327
## epsg (SRID): 26910
## proj4string: +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## name geometry
## 1 Space Needle POINT (548894.1 5274327)
# st_crs(hospitals)
# st_crs(zip_wa)
# st_crs(h2o)
h2o <- h2o %>% st_transform(4326)
hospitals <- hospitals %>% st_transform(4326)
Geoprocessing
Buffering
# create the points
snxy <- data.frame(name = "Space Needle", x = -122.3493, y = 47.6205)
space_needle <- st_as_sf(snxy, coords = c("x", "y"), crs = 4326)
shxy <- data.frame(name = "Savery Hall", x = -122.3083, y = 47.6572)
savery_hall <- st_as_sf(shxy, coords = c("x", "y"), crs = 4326)
zooxy <- data.frame(name = "Woodland Park Zoo", x = -122.3543, y = 47.6685)
wp_zoo <- st_as_sf(zooxy, coords = c("x", "y"), crs = 4326)
pts <- rbind(space_needle, savery_hall, wp_zoo)
# make the buffer with inline transforms
pts_buf_1km <- pts %>% st_transform(26910) %>% st_buffer(dist = 1000) %>% st_transform(2926)
# crs to be converted to Cartesian projected coordinate for the same measure of distance.
# write to the GPKG
mygpkg <- file.path(mydatadir, "r_gis.gpkg")
st_write(obj = pts_buf_1km, dsn = mygpkg, layer = "pts_buf_1km", quiet = TRUE, update = TRUE)
if(! exists("kctrans")){
kctrans <- st_read(
file.path(mydatadir, "Metro_Transportation_Network_TNET_in_King_County__trans_network_line.shp"),
quiet = TRUE)
}
# freeways are KC_FCC_ID = F
kcfwy <- kctrans %>% filter(KC_FCC_ID == "F")
# buffer
kcfwy_buf_500ft <- kcfwy %>% st_transform(2926) %>% st_buffer(500)
# write to the GPKG
mygpkg <- file.path(mydatadir, "r_gis.gpkg")
st_write(obj = kcfwy_buf_500ft, dsn = mygpkg, layer = "kcfwy_buf_500ft", quiet = TRUE, update = TRUE)
Point-in-polygon
# if no API key,
# acs5_2018_bg <- st_read(dsn = file.path(mydatadir, "census.gpkg"), layer = "acs5_2018_bg", quiet = TRUE)
# st_crs(acs5_2018_bg) <- 2926
# cache data
options(tigris_use_cache = TRUE)
# where to store data
tigris_cache_dir <- mydatadir
# if you have your API key, enter it here rather than using the system environment variable
# myapikey <- "foobar"
myapikey <- Sys.getenv("CENSUS_API_KEY")
census_api_key(myapikey)
## To install your API key for use in future sessions, run this function with `install = TRUE`.
# get the data and project it to match the bus stops, also calculate the area
acs5_2018_bg <- get_acs(
geography = "block group",
variables = c(medfamincome="B19113_001"),
state = "WA",
county = "King",
geometry = TRUE,
moe = 95,
cache_table = TRUE,
output = "wide") %>% # or tidy
st_transform(2926) %>%
mutate(area_ft = as.numeric(st_area(.)))
## Getting data from the 2013-2017 5-year ACS
colnames(acs5_2018_bg) <- tolower(colnames(acs5_2018_bg))
acs5_2018_bg %>%
ggplot() +
geom_sf(aes(fill = medfamincomee), size = .25) +
scale_fill_viridis_c() +
theme_void()

busstop <- st_read(
file.path(mydatadir, "busstop/busstop.shp"), quiet = TRUE)
st_crs(busstop) <- 2926
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that
colnames(busstop) <- tolower(colnames(busstop))
print(colnames(busstop))
## [1] "objectid" "busstop_id" "zonekey" "begin_date" "end_date"
## [6] "tlink_id" "ramkey" "access_lvl" "landing_pd" "trnsfr_pt"
## [11] "authority" "awning" "bearing_cd" "curb" "curb_ht"
## [16] "bay" "bike_rack" "created_by" "cross_stre" "curb_paint"
## [21] "dt_created" "dt_mapped" "dt_mod" "displacemt" "rampositio"
## [26] "ditch" "ext_surfc" "ext_width" "frm_cross" "frm_intrst"
## [31] "juris" "reg_fare" "rfp_dist" "zip_code" "i_sgn_anch"
## [36] "i_sgn" "int_loc" "link_len" "pcnt_from" "t_nd_from"
## [41] "mod_by" "news_box" "non_mt_sgn" "bollards" "num_shelt"
## [46] "on_street" "othr_cov_a" "owner" "paint_len" "pk_stp_sur"
## [51] "pullout" "ret_wall" "rfa_flag" "rt_sgn_tp" "sched_hold"
## [56] "shdr_surf" "shdr_width" "side" "side_cross" "side_on"
## [61] "swlk_width" "sgn_mt_dir" "sgn_pst_an" "sgn_pst_tp" "spc_sgn_tp"
## [66] "length" "status" "stop_type" "address" "add_commnt"
## [71] "strp_width" "t_signal" "wlk_surf" "xcoord" "ycoord"
## [76] "xcoord_off" "ycoord_off" "geometry"
busstop <- busstop %>% st_join(acs5_2018_bg)
print(colnames(busstop))
## [1] "objectid" "busstop_id" "zonekey" "begin_date"
## [5] "end_date" "tlink_id" "ramkey" "access_lvl"
## [9] "landing_pd" "trnsfr_pt" "authority" "awning"
## [13] "bearing_cd" "curb" "curb_ht" "bay"
## [17] "bike_rack" "created_by" "cross_stre" "curb_paint"
## [21] "dt_created" "dt_mapped" "dt_mod" "displacemt"
## [25] "rampositio" "ditch" "ext_surfc" "ext_width"
## [29] "frm_cross" "frm_intrst" "juris" "reg_fare"
## [33] "rfp_dist" "zip_code" "i_sgn_anch" "i_sgn"
## [37] "int_loc" "link_len" "pcnt_from" "t_nd_from"
## [41] "mod_by" "news_box" "non_mt_sgn" "bollards"
## [45] "num_shelt" "on_street" "othr_cov_a" "owner"
## [49] "paint_len" "pk_stp_sur" "pullout" "ret_wall"
## [53] "rfa_flag" "rt_sgn_tp" "sched_hold" "shdr_surf"
## [57] "shdr_width" "side" "side_cross" "side_on"
## [61] "swlk_width" "sgn_mt_dir" "sgn_pst_an" "sgn_pst_tp"
## [65] "spc_sgn_tp" "length" "status" "stop_type"
## [69] "address" "add_commnt" "strp_width" "t_signal"
## [73] "wlk_surf" "xcoord" "ycoord" "xcoord_off"
## [77] "ycoord_off" "geoid" "name" "medfamincomee"
## [81] "medfamincomem" "area_ft" "geometry"
# tabulate the count of transit stops
nbusstop <- busstop %>%
group_by(geoid) %>%
summarise(n_busstop = n(),
density_ha = n() / min(area_ft) * 107639 ,
medfamincomee = min(medfamincomee))
nbusstop %>% ggplot(aes(x = medfamincomee, y = density_ha)) +
geom_point() +
geom_smooth(method = "lm") +
xlab("block group median family income, ACS-5, 2018") + ylab("transit stop density per ha")
## Warning: Removed 66 rows containing non-finite values (stat_smooth).
## Warning: Removed 66 rows containing missing values (geom_point).

pander(
summary(
lm(data = nbusstop, medfamincomee ~ density_ha)))
| (Intercept) |
112170 |
1854 |
60.49 |
0 |
| density_ha |
89.16 |
10775 |
0.008275 |
0.9934 |
Fitting linear model: medfamincomee ~ density_ha
| 1198 |
47100 |
5.725e-08 |
-0.0008361 |
Polygon-on-polygon
nhood <- st_read(
file.path(mydatadir, "Community_Reporting_Areas.shp"),
quiet = TRUE) %>%
st_transform(2926)
names(nhood) = tolower(names(nhood))
# get the data and project it to match the bus stops, also calculate the area
acs5_2018_trt <- get_acs(
year = 2018,
geography = "tract",
variables = c(n="B06012_001", n_pov="B06012_002"),
state = "WA",
county = "King",
geometry = TRUE,
moe = 95,
cache_table = TRUE,
output = "wide") %>%
st_transform(2926) %>%
mutate(area_ft_tract = as.numeric(st_area(.)))
## Getting data from the 2014-2018 5-year ACS
colnames(acs5_2018_trt) <- tolower(colnames(acs5_2018_trt))
nhood_trt <- st_intersection(x = nhood, acs5_2018_trt) %>%
mutate(area_ft_intersect = as.numeric(st_area(.)),
n_est = ne * as.numeric(st_area(.)) / area_ft_tract,
n_est_pov = n_pove * as.numeric(st_area(.)) / area_ft_tract)
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries
nhood_pov <- nhood_trt %>%
group_by(gen_alias) %>%
summarize(
neighdist = first(neighdist),
n = sum(n_est),
n_pov = sum(n_est_pov),
pct_pov = round(sum(n_est_pov) / sum(n_est) * 100, 1))
nhood_pov %>%
ggplot() +
geom_sf(aes(fill = pct_pov), size = .25) +
scale_fill_viridis_c() +
theme_void()

nhood_pov %>%
ggplot(aes(x = reorder(gen_alias, pct_pov), y=pct_pov)) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
xlab("neighborhood") + ylab("% living under\nfederal poverty level")

st_write(obj = nhood_pov,
dsn = mygpkg,
layer = "nhood_pov",
quiet = TRUE,
update = TRUE,
delete_layer = TRUE)
st_write(obj =
st_union(kcfwy_buf_500ft),
dsn = mygpkg,
layer = "kcfwy_buf_500ft_union",
quiet = TRUE, update = TRUE)
# summarize == union
districts <- nhood_pov %>%
group_by(neighdist) %>%
summarise(
n = sum(n),
n_pov = sum(n_pov),
pct_pov = round(sum(n_pov) / sum(n) * 100, 1))
# save
st_write(obj = districts, dsn = mygpkg, layer = "districts", quiet = TRUE, update = TRUE, delete_layer = TRUE)
Leaflet map
mapview(nhood_pov, zcol = "pct_pov", legend = TRUE)
# CRS
nhood_pov_4326 <- nhood_pov %>% st_transform(4326)
# make a label with the count of persons, persons below poverty, and % poverty
nhood_pov_4326$mylab <- sprintf("n_pov=%s<br>n=%s<br>%s<br>",
round(nhood_pov_4326$n_pov, 0),
round(nhood_pov_4326$n, 0),
paste("pov=", nhood_pov_4326$pct_pov, "%", sep=""))
l <- leaflet(data = nhood_pov_4326) %>%
addPolygons(popup = ~mylab, weight = 2) %>%
addTiles()
l
# CRS
nhood_pov_4326 <- nhood_pov %>% st_transform(4326)
mypalette <- colorQuantile(palette = "viridis", domain = nhood_pov_4326$pct_pov, n = 4)
# make a label with the count of persons, persons below poverty, and % poverty
nhood_pov_4326$mylab <- sprintf("%s<br>%s<br>",
nhood_pov_4326$gen_alias,
paste("pov=", nhood_pov_4326$pct_pov, "%", sep=""))
l <- leaflet(data = nhood_pov_4326) %>%
addPolygons(popup = ~mylab,
weight = 2,
fillColor = ~mypalette(pct_pov),
opacity = 0.7) %>%
addTiles() %>%
addLegend(pal=mypalette, values=~pct_pov, opacity=0.7, title = "% below poverty", position = "bottomleft" )
l
# CRS
nhood_pov_4326 <- nhood_pov %>% st_transform(4326)
# hospitals
hospitals <- st_read(file.path(mydatadir, "medical_facilities", "medical_facilities.shp"), quiet = TRUE) %>% st_transform(4326)
mypalette <- colorQuantile(palette = "viridis", domain = nhood_pov_4326$pct_pov, n = 4)
# make a label with the count of persons, persons below poverty, and % poverty
nhood_pov_4326$mylab <- sprintf("%s<br>%s<br>",
nhood_pov_4326$gen_alias,
paste("pov=", nhood_pov_4326$pct_pov, "%", sep=""))
l <- leaflet(data = nhood_pov_4326) %>%
addPolygons(popup = ~mylab,
weight = 2,
fillColor = ~mypalette(pct_pov),
opacity = 0.7) %>%
addTiles() %>%
addLegend(pal=mypalette, values=~pct_pov, opacity=0.7, title = "% below poverty", position = "bottomleft" )
l <- addCircleMarkers(map = l,
data = hospitals,
radius = 5,
weight = 1,
opacity = 0.9,
fillOpacity = 0.5,
label = ~ABB_NAME)
l